# =============================================================================
# GHOST CELL MANAGEMENT AND LOW-LEVEL BOUNDARY IMPLEMENTATIONS
# =============================================================================

"""
    apply_periodic_bc!(field_data, grid, axis)

Apply periodic boundary conditions by copying opposite boundary values.
"""
function apply_periodic_bc!(field_data::AbstractArray, grid::AbstractGrid, axis::Int)
    # For periodic BC, copy values from opposite boundaries
    # This ensures continuity across the periodic domain
    
    dims = size(field_data)
    ndims_field = length(dims)
    
    if axis <= ndims_field
        # Get slice indices for boundaries
        lower_idx = [i == axis ? 1 : Colon() for i in 1:ndims_field]
        upper_idx = [i == axis ? dims[axis] : Colon() for i in 1:ndims_field]
        
        # Copy upper boundary to lower (enforce continuity)
        field_data[lower_idx...] = field_data[upper_idx...]
    end
    
    return field_data
end

"""
    apply_dirichlet_bc!(field_data, grid, axis, value)

Apply Dirichlet boundary condition with ghost cell update.
"""
function apply_dirichlet_bc!(field_data::AbstractArray, grid::AbstractGrid, axis::Int, value::Float64)
    # Apply Dirichlet BC to both boundaries
    apply_dirichlet_bc_side!(field_data, grid, axis, lower_side, value)
    apply_dirichlet_bc_side!(field_data, grid, axis, upper_side, value)
    return field_data
end

"""
    apply_dirichlet_bc_side!(field_data, grid, axis, side, value)

Apply Dirichlet BC to specific boundary side.
"""
function apply_dirichlet_bc_side!(field_data::AbstractArray, grid::AbstractGrid, axis::Int, side::BoundarySide, value::Float64)
    dims = size(field_data)
    ndims_field = length(dims)
    
    if axis <= ndims_field
        if side == lower_side
            # Set lower boundary to specified value
            idx = [i == axis ? 1 : Colon() for i in 1:ndims_field]
            field_data[idx...] .= value
        else  # upper_side
            # Set upper boundary to specified value
            idx = [i == axis ? dims[axis] : Colon() for i in 1:ndims_field]
            field_data[idx...] .= value
        end
    end
    
    return field_data
end

"""
    apply_neumann_bc!(field_data, grid, axis, derivative_value, order)

Apply Neumann boundary condition using ghost cells.
"""
function apply_neumann_bc!(field_data::AbstractArray, grid::AbstractGrid, axis::Int, derivative_value::Float64, order::Int=1)
    # Apply Neumann BC to both boundaries
    apply_neumann_bc_side!(field_data, grid, axis, lower_side, derivative_value, order)
    apply_neumann_bc_side!(field_data, grid, axis, upper_side, derivative_value, order)
    return field_data
end

"""
    apply_neumann_bc_side!(field_data, grid, axis, side, derivative_value, order)

Apply Neumann BC to specific boundary side using finite difference ghost cells.
"""
function apply_neumann_bc_side!(field_data::AbstractArray, grid::AbstractGrid, axis::Int, side::BoundarySide, derivative_value::Float64, order::Int=1)
    dims = size(field_data)
    ndims_field = length(dims)
    
    if axis <= ndims_field
        h = spacing(grid)[axis]
        
        if side == lower_side
            # Lower boundary: ∂f/∂n ≈ (f[2] - f[0]) / (2h) = derivative_value
            # Solve for ghost value: f[0] = f[2] - 2h * derivative_value
            interior_idx = [i == axis ? 2 : Colon() for i in 1:ndims_field]
            boundary_idx = [i == axis ? 1 : Colon() for i in 1:ndims_field]
            
            if dims[axis] >= 2
                # Use second-order finite difference
                field_data[boundary_idx...] = field_data[interior_idx...] .- 2 * h * derivative_value
            end
        else  # upper_side
            # Upper boundary: ∂f/∂n ≈ (f[n+1] - f[n-1]) / (2h) = derivative_value  
            # Solve for ghost value: f[n+1] = f[n-1] + 2h * derivative_value
            interior_idx = [i == axis ? dims[axis]-1 : Colon() for i in 1:ndims_field]
            boundary_idx = [i == axis ? dims[axis] : Colon() for i in 1:ndims_field]
            
            if dims[axis] >= 2
                # Use second-order finite difference
                field_data[boundary_idx...] = field_data[interior_idx...] .+ 2 * h * derivative_value
            end
        end
    end
    
    return field_data
end

"""
    apply_robin_bc!(field_data, grid, axis, alpha, beta, value)

Apply Robin boundary condition: α⋅f + β⋅∂f/∂n = value
"""
function apply_robin_bc!(field_data::AbstractArray, grid::AbstractGrid, axis::Int, alpha::Float64, beta::Float64, value::Float64)
    apply_robin_bc_side!(field_data, grid, axis, lower_side, alpha, beta, value)
    apply_robin_bc_side!(field_data, grid, axis, upper_side, alpha, beta, value)
    return field_data
end

"""
    apply_robin_bc_side!(field_data, grid, axis, side, alpha, beta, value)

Apply Robin BC to specific boundary side.
"""
function apply_robin_bc_side!(field_data::AbstractArray, grid::AbstractGrid, axis::Int, side::BoundarySide, alpha::Float64, beta::Float64, value::Float64)
    dims = size(field_data)
    ndims_field = length(dims)
    
    if axis <= ndims_field && dims[axis] >= 2
        h = spacing(grid)[axis]
        
        if side == lower_side
            # Robin BC at lower boundary: α⋅f[1] + β⋅(f[2] - f[0])/(2h) = value
            # Solve for ghost value f[0]: f[0] = f[2] - (2h/β)(value - α⋅f[1])
            if abs(beta) > 1e-12
                boundary_idx = [i == axis ? 1 : Colon() for i in 1:ndims_field]
                interior_idx = [i == axis ? 2 : Colon() for i in 1:ndims_field]
                ghost_idx = boundary_idx  # Use boundary as ghost approximation
                
                f_boundary = field_data[boundary_idx...]
                f_interior = field_data[interior_idx...]
                
                # Ghost value that satisfies Robin BC
                field_data[ghost_idx...] = f_interior .- (2 * h / beta) .* (value .- alpha .* f_boundary)
            else
                # Pure Dirichlet case (β=0)
                apply_dirichlet_bc_side!(field_data, grid, axis, side, value / alpha)
            end
        else  # upper_side
            # Similar for upper boundary
            if abs(beta) > 1e-12
                boundary_idx = [i == axis ? dims[axis] : Colon() for i in 1:ndims_field]
                interior_idx = [i == axis ? dims[axis]-1 : Colon() for i in 1:ndims_field]
                ghost_idx = boundary_idx
                
                f_boundary = field_data[boundary_idx...]
                f_interior = field_data[interior_idx...]
                
                field_data[ghost_idx...] = f_interior .+ (2 * h / beta) .* (value .- alpha .* f_boundary)
            else
                apply_dirichlet_bc_side!(field_data, grid, axis, side, value / alpha)
            end
        end
    end
    
    return field_data
end

"""
    apply_bc_generic!(field_data, bc, grid, axis, time)

Generic fallback for boundary condition application.
"""
function apply_bc_generic!(field_data::AbstractArray, bc::BoundaryCondition, grid::AbstractGrid, axis::Int, time::Float64)
    @warn "Using generic BC application for $(typeof(bc)) - may not be implemented"
    return field_data
end

# =============================================================================
# BOUNDARY CONDITION PARSING AND STRING INTERFACE
# =============================================================================

"""
    parse_boundary_condition(bc_string::String) -> BoundaryCondition

Parse boundary condition from string representation.

# Supported formats:
- "periodic" -> PeriodicBC()
- "dirichlet:value" -> DirichletBC(value)  
- "neumann:value" -> NeumannBC(value)
- "robin:alpha,beta,value" -> RobinBC(alpha, beta, value)

# Examples:
```julia
bc = parse_boundary_condition("dirichlet:1.0")
bc = parse_boundary_condition("neumann:0.0")  
bc = parse_boundary_condition("robin:1.0,0.1,0.0")
```
"""
function parse_boundary_condition(bc_string::String)::BoundaryCondition
    bc_string = strip(lowercase(bc_string))
    
    if bc_string == "periodic"
        return PeriodicBC()
    elseif startswith(bc_string, "dirichlet:")
        value_str = bc_string[11:end]  # Remove "dirichlet:"
        value = parse(Float64, value_str)
        return DirichletBC(value)
    elseif startswith(bc_string, "neumann:")
        value_str = bc_string[9:end]  # Remove "neumann:"
        value = parse(Float64, value_str)
        return NeumannBC(value)
    elseif startswith(bc_string, "robin:")
        params_str = bc_string[7:end]  # Remove "robin:"
        params = split(params_str, ",")
        length(params) >= 2 || throw(ArgumentError("Robin BC requires at least alpha,beta parameters"))
        
        alpha = parse(Float64, strip(params[1]))
        beta = parse(Float64, strip(params[2]))
        value = length(params) >= 3 ? parse(Float64, strip(params[3])) : 0.0
        
        return RobinBC(alpha, beta, value)
    else
        throw(ArgumentError("Unknown boundary condition format: $bc_string"))
    end
end

"""
    parse_boundary_conditions(bc_strings::Vector{String}) -> Vector{BoundaryCondition}

Parse multiple boundary conditions from string vector.
"""
function parse_boundary_conditions(bc_strings::Vector{String})::Vector{BoundaryCondition}
    return [parse_boundary_condition(bc_str) for bc_str in bc_strings]
end

# =============================================================================
# GHOST CELL MANAGEMENT AND CACHING
# =============================================================================

"""
    get_ghost_values(bc::BoundaryCondition, field_data, grid, axis, side)

Get ghost cell values for boundary condition.
"""
function get_ghost_values(bc::BoundaryCondition, field_data::AbstractArray, grid::AbstractGrid, axis::Int, side::BoundarySide)
    # Default implementation - copy boundary values
    dims = size(field_data)
    
    if side == lower_side
        idx = [i == axis ? 1 : Colon() for i in 1:length(dims)]
    else
        idx = [i == axis ? dims[axis] : Colon() for i in 1:length(dims)]
    end
    
    return field_data[idx...]
end

function get_ghost_values(bc::DirichletBC, field_data::AbstractArray, grid::AbstractGrid, axis::Int, side::BoundarySide)
    # For Dirichlet BC, ghost values enforce the boundary value
    value = evaluate_bc_value(bc.value, grid, axis, 0.0)
    dims = size(field_data)
    
    # Create array of appropriate size filled with boundary value
    if side == lower_side
        idx = [i == axis ? 1 : Colon() for i in 1:length(dims)]
    else
        idx = [i == axis ? dims[axis] : Colon() for i in 1:length(dims)]
    end
    
    ghost_shape = [i == axis ? 1 : dims[i] for i in 1:length(dims)]
    return fill(value, ghost_shape...)
end

function get_ghost_values(bc::NeumannBC, field_data::AbstractArray, grid::AbstractGrid, axis::Int, side::BoundarySide)
    # For Neumann BC, ghost values extrapolate to maintain desired derivative
    value = evaluate_bc_value(bc.value, grid, axis, 0.0)
    h = spacing(grid)[axis]
    dims = size(field_data)
    
    if side == lower_side && dims[axis] >= 2
        interior_idx = [i == axis ? 2 : Colon() for i in 1:length(dims)]
        return field_data[interior_idx...] .- 2 * h * value
    elseif side == upper_side && dims[axis] >= 2
        interior_idx = [i == axis ? dims[axis]-1 : Colon() for i in 1:length(dims)]
        return field_data[interior_idx...] .+ 2 * h * value
    else
        # Fallback
        boundary_idx = [i == axis ? (side == lower_side ? 1 : dims[axis]) : Colon() for i in 1:length(dims)]
        return field_data[boundary_idx...]
    end
end

# =============================================================================
# INTEGRATION WITH DIFFERENTIAL OPERATORS
# =============================================================================

"""
    modify_operator_for_bc!(op_matrix, bc::BoundaryCondition, grid, axis)

Modify differential operator matrix to incorporate boundary conditions.
"""
function modify_operator_for_bc!(op_matrix::AbstractMatrix, bc::BoundaryCondition, grid::AbstractGrid, axis::Int)
    # Default: no modification
    return op_matrix
end

function modify_operator_for_bc!(op_matrix::AbstractMatrix, bc::DirichletBC, grid::AbstractGrid, axis::Int)
    # For Dirichlet BC, modify operator to enforce boundary values
    n = size(op_matrix, 1)
    dims = size(grid)
    
    if axis <= length(dims)
        # Set boundary rows to identity (for homogeneous BC)
        # This is a simplified version - full implementation would be more complex
        op_matrix[1, :] .= 0.0
        op_matrix[1, 1] = 1.0
        
        if n > 1
            op_matrix[n, :] .= 0.0
            op_matrix[n, n] = 1.0
        end
    end
    
    return op_matrix
end

function modify_operator_for_bc!(op_matrix::AbstractMatrix, bc::PeriodicBC, grid::AbstractGrid, axis::Int)
    # For periodic BC, connect boundaries in operator
    n = size(op_matrix, 1)
    dims = size(grid)
    
    if axis <= length(dims) && n >= 2
        # Connect first and last rows to enforce periodicity
        # This is simplified - full implementation would handle multi-dimensional indexing
        op_matrix[1, n] += op_matrix[1, 1]
        op_matrix[n, 1] += op_matrix[n, n]
    end
    
    return op_matrix
end

# =============================================================================
# UTILITY FUNCTIONS
# =============================================================================

"""
    is_homogeneous(bc::BoundaryCondition) -> Bool

Check if boundary condition is homogeneous (zero value).
"""
function is_homogeneous(bc::BoundaryCondition)
    return false  # Default: not homogeneous
end

function is_homogeneous(bc::DirichletBC)
    return bc.homogeneous
end

function is_homogeneous(bc::NeumannBC)
    return isa(bc.value, Float64) && bc.value ≈ 0.0
end

function is_homogeneous(bc::RobinBC)
    return isa(bc.value, Float64) && bc.value ≈ 0.0
end

"""
    bc_string(bc::BoundaryCondition) -> String

Convert boundary condition to string representation.
"""
function bc_string(bc::PeriodicBC)
    return "periodic"
end

function bc_string(bc::DirichletBC)
    if isa(bc.value, Float64)
        return "dirichlet:$(bc.value)"
    else
        return "dirichlet:<function>"
    end
end

function bc_string(bc::NeumannBC)
    if isa(bc.value, Float64)
        return "neumann:$(bc.value)"
    else
        return "neumann:<function>"
    end
end

function bc_string(bc::RobinBC)
    if isa(bc.value, Float64)
        return "robin:$(bc.alpha),$(bc.beta),$(bc.value)"
    else
        return "robin:$(bc.alpha),$(bc.beta),<function>"
    end
end

function bc_string(bc::MixedBC)
    return "mixed:[$(bc_string(bc.lower)), $(bc_string(bc.upper))]"
end

"""
    summarize_boundary_conditions(grid::AbstractGrid)

Print summary of boundary conditions for a grid.
"""
function summarize_boundary_conditions(grid::AbstractGrid)
    bcs = get_boundary_conditions(grid)
    println("Boundary conditions for $(typeof(grid).name.name):")
    for (i, bc) in enumerate(bcs)
        println("  Dimension $i: $(bc_string(bc))")
    end
end